import numpy as np
import struct
import timeit

## 非压缩格式得带状矩阵求解
def gaussin1(a,b,n,p,q):
    ## 系数矩阵a，右端向量b，n阶矩阵，上带宽q,下带宽p
    print("正在消去...")
    cout = 0  # 定义计算次数
    for k in range(n - 1):  # k表示第一层循环，(1，n-1)行
        # 限制条件
        if((k+p) < n):
            mm = k+p+1
        else:
            mm = n
        if((k+q) < n):
            nn = k+q+1
        else:
            nn = n
        if (a[k][k] == 0):
            print("no answer")
        for i in range(k + 1, mm):  # i表示第二层循环,(k+1,mm)行,计算该行消元的系数
            l = a[i][k] / a[k][k]
            cout += 1
            for j in range(k, nn):  # 第i行进行消去
                a[i][j] -= l * a[k][j]
                cout += 1
            b[i] -= l * b[k]
            cout += 1
    # 回代求出方程解
    print("回代求出方程解...")
    print("显示前20位x[i]...")
    x = np.zeros(n)
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]  # 先算最后一位的x解
    for i in range(n - 2, -1, -1):  # 依次回代倒着算每一个解
        if((i+q)<n):
            nn1 = i+q+1
        else:
            nn1 = n
        for j in range(i + 1, nn1):
            b[i] -= a[i][j] * x[j]
            cout += 1
        x[i] = b[i] / a[i][i]
        cout += 1
    for i in range(20):  # 显示前20位计算解，保留小数点后三位
        print("x[" + str(i + 1) + "] = ", np.around(x[i],3))
    print("浮点数计算次数约为：", cout)
    return 0

## 压缩格式得带状矩阵求解
def gaussin2(a,b,n,p,q):
    ## 系数矩阵a，右端向量b，n阶矩阵，上带宽q,下带宽p
    cout = 0
    print("正在消去...")
    for k in range(n-1):
        if ((k+p)<n):
            mm = k+p+1
        else:
            mm = n
        nn = p-1
        for i in range(k+1,mm):
            l = a[i][nn]/a[k][p]
            cout += 1
            mn = p + 1
            for j in range(nn+1, nn+q+1):
                a[i][j] -= l*a[k][mn]
                cout += 1
                mn += 1
            b[i] -= l*b[k]
            cout += 1
            nn -= 1
    # 回代求出方程解
    print("回代求出方程解...")
    print("显示前20位x[i]...")
    x = np.zeros(n)
    x[n - 1] = b[n - 1] / a[n - 1][p]  # 先算最后一位的x解
    for i in range(n - 2, -1, -1):  # 依次回代倒着算每一个解
        m1 = i+1
        for j in range(p+1, p+q+1):
            if (m1 < n):
                b[i] -= a[i][j] * x[m1]
                cout += 1
            m1 += 1
        x[i] = b[i] / a[i][p]
        cout += 1
    for i in range(20):  # 显示前20位计算解，保留小数点后三位
        print("x[" + str(i + 1) + "] = ", np.around(x[i], 3))
    print("浮点数计算次数约为：", cout)
    return 0
## 读文件，开始解方程
def readfile(FName):
    start = timeit.default_timer()
    typename = {0x102: '非压缩格式', 0x202: '压缩格式'}
    print("\n\n数据文件:", FName)
    f = open(FName, 'rb')
    id, ver, id1 = struct.unpack('<3I', f.read(12))  # 读入文件头标志部分
    n, q, p = struct.unpack('<3I', f.read(12))  # 读入系数矩阵结构信息
    print('id:', hex(id), '  ver:', hex(ver), '  id1:', hex(id1), "   ", typename[ver])
    print('阶数:', n, ' 上带宽:', q, ' 下带宽:', p)

    B = np.zeros([n])  # 方程右端向量
    if hex(ver) == '0x202':  # 压缩格式求解带状矩阵
        m = p + q + 1
        A = np.zeros([n, m])
        print("读入系数矩阵...")
        for i in range(n):
            for j in range(m):
                A[i, j], = struct.unpack('f', f.read(4))
        print("读入右端常量...")
        for i in range(n):
            B[i], = struct.unpack('f', f.read(4))
        gaussin2(A, B, n, p, q)
    else:  # 非压缩格式求解带状矩阵
        A = np.zeros([n, n])
        print("读入系数矩阵...")
        for i in range(n):
            for j in range(n):
                A[i, j], = struct.unpack('f', f.read(4))
        print("读入右端常量...")
        for i in range(n):
            B[i], = struct.unpack('f', f.read(4))
        gaussin1(A, B, n, p, q)
    f.close()
    end = timeit.default_timer()
    tmp = np.around(1000 * float(end - start), 2)
    print("计算时长 =", tmp, "ms")

if __name__ == '__main__':
    # FName = "C:/Users/yc_scorpio/Desktop/test/data20201.dat"
    # FName = "C:/Users/yc_scorpio/Desktop/test/data20202.dat"
    # FName = "C:/Users/yc_scorpio/Desktop/test/data20203.dat"
    # FName = "C:/Users/yc_scorpio/Desktop/test/data20204.dat"
    FName = "C:/Users/yc_scorpio/Desktop/test/data20205.dat"
    readfile(FName)